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We provide strong evidence that, up to 99.999% of extremality, Kerr-Newman black holes (KN 
BHs) are linear mode stable within Einstein-Maxwell theory. We derive and solve, numerically, a 
coupled system of two PDEs for two gauge invariant fields that describe the most general linear 
perturbations of a KN BH (except for trivial modes that shift the parameters of the solution). 
We determine the qnasinormal mode (QNM) spectrum of the KN BH as a function of its three 
parameters and find no unstable modes. In addition, we find that the QNMs that are connected 
continuously to the gravitational i = m = 2 Schwarzschild QNM dominate the spectrum for all 
values of the parameter space (m is the azimuthal number of the wave function and £ measures 
the number of nodes along the polar direction). Furthermore, all QNMs with i = m approach 
Reoj = and Imo; = 0 at extremality; this is a universal property for any field of arbitrary 

spin |s| < 2 propagating on a KN BH background {uj is the wave frequency and the BH angular 
velocity at extremality). We compare our results with available perturbative results in the small 
charge or small rotation regimes and find good agreement. We also present a simple proof that the 
Regge-Wheeler (odd) and Zerilli (even) sectors of Schwarzschild perturbations must be isospectral. 


Introduction. The uniqueness theorems mm state 
that the Kerr-Newman black hole (KN BH) [31 Sj is 
the unique, most general family of stationary asymptot¬ 
ically flat BHs, of Einstein-Maxwell theory. It is char¬ 
acterised by 3 parameters: mass M, angular momentum 
J = Ma and charge Q. The Kerr, Reissner-Nordstrom 
(RN) and Schwarzschild (Schw) BHs constitute limiting 
cases: Q = 0, a = 0 and Q = a = 0, respectively. 

Given their uniqueness, the most relevant question to 
consider is the linear mode stability of these BHs. It 
is known that the Kerr, RN, and Schw BHs are linear 
mode stable. Indeed, the perturbation study of the lin¬ 
earised Einstein(-Maxwell) equation gives the quasinor¬ 
mal mode (QNM) spectrum of frequencies (that describes 
the damped oscillations of the BH back to equilibrium) 
resulting in no unstable modes [5H2D] (see review [2T]). 
Remarkably, for these BHs, the QNM spectrum turns 
out be encoded in a single separable equation—known 
as the Regge-Wheeler-Zerilli equation [5H7] (for RN and 
Schw) and the Teukolsky equation [T3] (for Kerr, RN and 
Schw)—that effectively yields a pair of ODEs. 

Unfortunately, it does not seem possible to cast a gen¬ 
eral perturbation of a KN BH as a single PDE. Therefore, 
obtaining the QNM spectrum of KN BHs requires solv¬ 
ing coupled PDEs. Naively, one expects to find a system 
of nine coupled PDEs. However, working in the so-called 
phantom gauge, = 0, Ghandrasekhar, re¬ 

duced the problem to the study of two coupled PDEs [13] 
(see also [52|). Despite this significant progress, finding 
the QNM spectrum and addressing the problem of the 
linear mode stability of the KN BH has remained a ma¬ 
jor open problem of Einstein-Maxwell theory since the 


80’s, when Chandrasekhar’s seminal work m was pub¬ 
lished. 

Recently there have been some notable efforts to ad¬ 
dress this problem. Refs. [231 El] and [22] have found 
the QNM spectrum in a perturbative small rotation and 
charge, respectively, expansion around the RN and Kerr 
BHs. These works find no sign of linear instability; how¬ 
ever such an instability is more likely to be found in ex¬ 
treme regimes where both Q and a are large. Another 
remarkable effort to infer the (non-)linear stability of KN 
BHs has been made in [35], where the full time evolu¬ 
tion of some KN BH with a given initial perturbation 
is considered, finding no sign of a nonlinear instability. 
However, since non-linear simulations are computation¬ 
ally costly, the search in moduli space is modest. 

In this letter, we derive two coupled PDEs that reduce 
to the Chandrasekhar coupled PDE system upon gauge 
fixing and compute the QNM spectrum of the KN BH 
to a high degree of accuracy. Up to ajaext = 0.99999 
we find no sign of a linear mode instability for any of 
the gravito-electromagnetic modes that are described by 
i = 1,2,3,4 and \m\ < t. We use two distinct numeri¬ 
cal methods that have been developed to solve efficiently 
similar systems of (several coupled) ODEs and PDEs that 
appear in QNM, superradiant and ultraspinning instabil¬ 
ity studies [35H55] . One of these methods formulates the 
problem as a quadratic eigenvalue problem in the fre¬ 
quency and employs a pseudospectral grid collocation. 
The other method searches directly for specific QNMs 
using a Newton-Raphson root-finding algorithm. We re¬ 
fer the reader to [26ll35] for details. The pseudospectral 
exponential convergence of our method, and the use of 


2 


quadruple precision, guarantees that the results are ac¬ 
curate up to, at least, one tenth of a decimal place. 

Notation: We use the standard nomenclature of the 
Newman-Penrose formalism to denote components of the 
curvatures, electromagnetic field strength and connec¬ 
tions [32]. denotes background quantities, while 

denotes a perturbed quantity at the linear order. 

Formulation of the problem. We write the KN 
BH solution in standard Boyer-Lindquist coordinates 
{t, r, 0, (f)} (time, radial, polar, azimuthal coordinates) 
[3|. Its event horizon, with angular velocity fin and 
temperature Tf/, is generated by the Killing vector K = 
dt + flndcj). The location of the horizon r_|_ is given by 
the largest root of the function A. These quantities are 
given in terms of the parameters {M, a, Q} as follows: 


A = r^ -2Mr + a^+ Q^, 


Hh = 


Th = 


r+ = M + -y/ 

1 

47rr_|_ -|- 
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The KN BH has a regular extremal configuration when 
its temperature vanishes and its angular velocity reaches 
a maximum. For fixed M and Q this occurs for a = 
ttext = Thus, at extremality we have 


with f = r + fa cos 0. We emphasise that 'ijj -2 and ip-i 
(as well as (p -2 and (p-i) are gauge invariant perturbed 
quantities, i.e. they are invariant under both linear dif- 
feomorphisms and tetrad rotations. Furthermore, these 
are the NP scalars that are relevant for the study of per¬ 
turbations that are outgoing at future null infinity and 
regular at the future horizon |44j . Fixing a gauge in which 
= 0, we obtain the Chandrasekhar coupled 
PDF system [T5| (see also the derivation in [22]). Fi¬ 
nally, note that in the limit Q 0 and/or a —)■ 0 these 
equations decouple yielding the Teukolsky equation. 

In order to solve these equations, we need to impose 
appropriate boundary conditions (BCs). The t — <j) sym¬ 
metry of the KN BH guarantees that we can consider 
only modes with m > 0, say, as long as we consider both 
positive and negative Re(a;); when a = 0 this enhances 
to a t —>■ —t symmetry and the QNM frequencies form 
pairs of {oj, —w*}. 

At spatial infinity, a Frobenius analysis of Q and the 
requirement that we have only outgoing waves fixes the 
decay to be (s = —2, —1) 
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We consider the most general perturbation of a KN 
BH. Using the fact that dt and are Killing vector fields 
of the KN background, we Fourier decompose its per¬ 
turbations as This introduces the frequency 

w and azimuthal quantum number m of the perturba¬ 
tion. By formulating the perturbation problem in the 
Newman-Penrose (NP) formalism, we obtain a set of two 
coupled partial differential equations that describe the 
most general linear perturbation of a KN BH (see Ap¬ 
pendix for details of the derivation) 

(0_2 + ^f^r- 2 ) <P-2 + Q-2V-1 = 0 , (3) 

(o_i + = 0, 

where differential operators {0,V,Q} are given in ^ of 
Appendix, ip -2 = and ip-i = 

Substituting the background values of the NP quanti¬ 
ties, the above equations reduce to 


{j'-2 + Q^G- 2 ) ^-2 + Q‘^'H-2'^-1 = 0 , (4) 

(j-_i + Q^g.i) V'-i + = 0, 

where second order differential operators {J^, 5,'H} are 
given in Q of Appendix and 

V'-2 = 


where as(d) and Ps{0) are related by two Robin BCs 
which are fixed by Q. 

At the horizon, a Frobenius analysis and requiring only 
regular modes in the ingoing Eddington-Finkelstein co¬ 
ordinates yields the near-horizon expansion 

- r+)“® [as(6») +6s(6>)(r-r+) H-], 

where bs{6) is related to as{0) and its derivatives. 

At the North (South) pole x = cos 9 = 1 (—1), regular¬ 
ity dictates that the fields behave as (e = 1 for m > 2, 
while e = — 1 for m = 0,1 modes) 

where Aj(r) and Bf{r) are related by two Robin BCs 
that are fixed by the equations of motion. 

We consider only modes with the lowest radial over¬ 
tone {n = 0) because these are the ones that have smaller 
|Ima;| and thus they are the ones that dominate a time 
evolution and are more likely to become unstable near 
extremality. Note also that we can scale out one of the 3 
parameters of the solution. Thus, we work with the adi- 
mensional parameters {a/M,Q/M} (or {a/r+,Q/r+}) 
and ujM. 

Results and Discussion. Our primary aim is to find 
whether KN BHs can be linear mode unstable. For that, 
we study the QNM spectrum and check if there are modes 
with Imw > 0. Note that for Q,a ^ 0 we ought to re¬ 
cover the Schw QNMs. In this limit, there are two fam¬ 
ilies of QNMs, namely the Regge-Wheeler (odd or ax¬ 
ial) modes and the Zerilli (even or polar) modes. These 
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FIG. 1: All lowest radial overtone QNMs of Q = a KN BHs FIG. 2: Real part of the QNM frequencies with I! = m = 2 of 
that start at the £ = 3 Schw gravitational QNM (red disc). KN BHs with Q = a and fixed Qlr+ = 0.0,0.2, • ■ • , 0.9 that 

start at the Schw gravitational QNM with £ = 2 (red disc). 


families are isospectral, i.e. they have exactly the same 
spectrum m-. more later). Thus we only need to dis¬ 
tinguish the gravitational modes (described in Table V 
of page 262 m —hereafter Table of m —by the eigen¬ 
function Z 2 ) from the electromagnetic modes (described 
in Table of [H] by the eigenfunction Zi). Each of these 
is specified by the harmonic number £ = 1, 2, 3, • • • {Z 2 
modes with I = \ are pure gauge modes). When the 
BH has charge and rotation, we have to scan a two pa¬ 
rameter space in {Q/M,a/M}. The above two families 
become coupled gravito-electromagnetic QNMs and the 
Schwarzschild eigenvalue £ does not appear explicitly in 
the KN PDEs Q. However, we can still count the num¬ 
ber of nodes along the polar direction of the eigenfunc¬ 
tions of Q and this gives £. 

We perform a complete scan in {Q/M,a/M} for all 
modes with £ = 1,2,3, \m\ < 3 (both in the Zi and Z 2 
sectors). Modes with £ = 4 are also studied, but there we 
focus on modes that approached Im a; = 0 at extremality. 
As one of our main results, we do not find any unstable 
mode with Imw > 0, even when we probe regions in 
parameter space for which ajoext = 0.99999. We see this 
as good numerical evidence that the KN BH is linearly 
mode stable. 

To illustrate our search, in Fig. [l]we take KN BHs with 
Q = a and we display all the QNMs that are continuously 
connected to the gravitational Z 2 Schw QNM with £ = 
3, namely uiM = 0.59944329 — 0.09270305 f (see Table 
of [H]). The different QNMs are distinguished by their 
azimuthal number m = 0,1,2, 3 and by whether they 
have positive or negative Re w (modes with m = 0 have a 
pair of QNMs {w, —w*}; see discussion above). All these 
modes become degenerate in the Schwarzschild limit (red 
disk in Fig. [^. We plot the imaginary (main plot) and 
real (inset plot) parts of the frequency ujM as a function 
of a/oext = aj^/M"^ — Q^. We see that the most likely 
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FIG. 3: Similar to Fig. [^but now for Im(tjM). 


mode to be unstable is the £ = 3 mode with Rew > 0 
(magenta triangles). However, we follow this mode up to 
a/oext = 0.99999 and find that although the Imw quickly 
approaches zero as a —>■ Uext, it never crosses Imw = 0. 

It is also relevant to ask what are the dominant QNMs, 
i.e. the modes with the smallest |ImtLi|. We find 
that the QNM family that, in the Q,a ^ 0 limit, ap¬ 
proaches the Z 2 Schw QNM with £ = m = 2 with 
ioM = 0.37367168 - 0.08896232f (Table of [H]) is the 
one that always (i.e. for a given Q and a) has the small¬ 
est |Ima;|. Therefore, these QNMs must be the dominant 
modes in a time evolution of the KN BH. Since this mode 
is the most relevant in a time evolution process, hereafter 
we will use it to illustrate our discussions (other modes 
will be presented elsewhere). 

The plots of Figs. (real part of wM) and[^ (imaginary 
part of ujM) give details of the Z 2 ,£ = m = 2 mode. 
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We represent the QNMs of the Q = a KN BH by blue 
disks, but we also present the QNMs for KN BHs with 
fixed charge Q/r^ (see plot legends) as the rotation grows 
from zero to aext/M. Fig. shows that as extremality is 
approached we always have Rew —>■ 20“*. On the other 
hand, Fig. shows that Imcj —0“ as extremality is 
approached . Again, we emphasise that the last point 
of each of these curves is, at least, at a = 0.99999 Oeait. 
Figs. Bi and 1^ illustrate a general property of the KN 
QNMs with i = m (both in the Z-^ and sectors): as 
extremality is approached one has Rew —)■ and 

Imw —)■ 0“. Collecting previous results m uni 133 EH, 
we can can now state that this is a universal property 
for any perturbation spin s. As we discussed above, in a 
3-dimensional {Q/M,a/M,lm{ujM)} plot, all these £ = 
m ^ 2 QNM are below the Z 2 , £ = m = 2. 

Previously, there were some attempts to find the 
QNMs of the KN BH using a perturbative analysis, no¬ 
tably for small a around the RN QNMs [231 Elj and for 
small Q around the Kerr QNMs [22] • We use these per¬ 
turbative results to further check our results for small a or 
Q, thus establishing the regimes of validity of the afore¬ 
mentioned approximations. To compare the two pertur¬ 
bative analyses in a single graphic it is convenient to look 
at QNMs with Q = a, see Figs. |^(real part) andjsUimag- 
inary part). We see that the approximations of[23l l24] 
([22] ) are within 1% of the exact results up to ~ 25% 
(~ 70%) of extremality. 



FIG. 4: Comparison (for Reoj) between the exact i = m = 2 
Z 2 QNMs of KN with Q = a (blue disks) with the small a 
approximations of [23] (red diamonds with their 1% error bar) 
and with the small Q approximations of [22] (green circles). 

Ref. [2H considered the full time evolution of some 
KN BHs with a given initial perturbation and found no 
sign of a nonlinear instability, which is consistent with 
our full parameter scan of the QNMs up to ajaext = 
0.99999. Ref. [23 also finds numerical evidence that 



FIG. 5: Similar comparison as in Fig. |^but now for Imoj. 

some £ = m = 2 QNMs of a KN BH (with a/Q > 1) 
should have the scaling uj = lo {^a/^J m'^ — . We can 

test this claim with very high accuracy and we find that 
it does not hold (although it must be emphasised that 
our linear results are well within the numerical accuracy 
of [23; they differ by at most 1%). Indeed, in Fig. 
we pick some KN BHs with fixed Q/r+ and a/Q > 1. 
For £ = m = 2, we plot the difference A{ujM) between 
ooM and the QNM frequency of the <5 = 0 KN BH as a 
function of a/ (real part of uiM in the main 
plot and imaginary part in the inset plot). If the scaling 
proposed in [23133 were present, all these curves should 
superpose at A{ujM) = 0. This is not the case: e.g. 
|AIm(a;M)| can be as high as 0.02 (when the error in our 
results is ±10“^°). Coming back to the fact that these 
modes approach at extremality, this had to be 

the case since {aj ( see extremal 

points on Fig. [^. To sum, although the proposed scaling 
fails to hold just by a small relative amount (less than 
1%), our numerical and analytical considerations show 
that it is only approximate, but not exact. 

Finally, we present a simple proof of the isospectral 
property m of the Schw and RN QNMs [33 ■ We can use 
the Teukolsky equation to study the QNMs of the Schw 
BH, instead of the Regge-Wheeler—Zerilli (RWZ) formal¬ 
ism. The two must give the same spectrum. The Teukol¬ 
sky formulation has a single gauge invariant variable ip -2 
that must translate into two gauge invariant variables in 
the RWZ formulation, namely Regge-Wheeler’s 4>v and 
Zerilli’s <i)s eigenfunctions. Refs. [33 131] 03 
unique differential map that applied to 'ip -2 constructs 
$v and another that constructs $s. Isospectrality is the 
statement that $v and <i>s have the same QNM spectrum. 
Since <i)v and <i)s are constructed from the same Teukol¬ 
sky ip- 2 , it follows that this must be the case. 
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■e- Qfr,=0.0 
Qfr,=0.1 
Q/r,=0.2 
-»- Q/r,=0.3 
Qlr^=O.A 
Q/r,=0.5 
* Q/r,=0.6 


ant quantity that comes out of the first equation, i.e. we 
manipulate the NP equations in such a way as to obtain 
an expression in which the only second order in derivative 
terms are of the form {DA — Subsequently, we 

simplify the resulting expression as before. The fact that 
the Maxwell NP eqns. (7.23) and (7.25) mix background 
and perturbed quantities does not cause any problems as 
the background contributions cancel once the expression 
is fully simplified as described above. 

Finally, the differential operators {0,V, Q} introduced 
in ^ are: 


FIG. 6: Some of the £ = m = 2 QNMs of Figs. and 
but truncated for a/Q > 1. They do not display the scaling 
conjectured in |25| . 
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0-2 = {A + 3-f--f* +ifi + ti*){D-p) 

-{6* -b 3 q! + -b 47r - T*){S -b 4/3 - r) - 34-2, 
r-2 = 2-A{A + 3'y-'y*)A-{D- p) 

—4(r* — 7r)A+(d* -b 4/3 — r), 

Q-2 = 2/$^°^*{(A-b37-7*-b2^)A_(r-b2a-b67r) 
-b(T* - 7r)A+(A -b 27 -b 6/r)}, 

0-1 = {A + Sj + j*+ 5p + p*){D-4p) (6) 

-((5* +a + /3* + 5n- t*){S -b 2/3 - 4 t), 

V-i = 2{D-Ap + p*)A+{A + 2-f + 6p) 

+2{d -b 3/3 - a* - 4r - 7r*)A_ (d* + 2a -b 67r), 
Q_i = -4$®{(3^-2p + p*)A+(r+4/3-T) 

+ ((5 -b 3/3 — a* — 2 t -b 7r*)A_(A -b 47 — p)}. 


where A± = (34'2°^ =F 2$^°^) 


Supplementary material: Derivation of the coupled 
equations Q 

The derivation of the two coupled PDEs is simi¬ 
lar to the derivation of the Chandrasekhar coupled PDEs 
[in]([i2]), except that we do not fix gauge. We derive the 
equations using the NP formalism. All the NP equations 
quoted in this section refer to equations in [36] . The first 
equation is derived by considering (5*(7.32d) — A(7.32c). 
The resulting expression is simplified by using these as 
well as other NP equations to express all perturbed quan¬ 
tities in terms of and In particular, NP 

eqn. (7.21j) is useful. The derivation of the second equa¬ 
tion is slightly less straightforward. Here we consider 
an appropriate combination of [73(7.32d) — d(7.32c)] and 
[A(7.23) —(5*(7.25)] that is suggested by the gauge invari- 


Supplementary material: The second order 
differential operators {£F, Q, 

The coupled system of two PDEs Q that describe the 
most general linear perturbation of a KN BH (in the 
Newman-Penrose formalism) were first derived, in the 
phantom gauge = 0, by Ghandrasekhar in 

his celebrated textbook |15j . A nice compact summary 
of this derivation can also be found in Appendix of |22j . 
This still yields our equations Q but this time with the 
associated gauge fixing applied to the variables (§. 

The second order differential operators {J^, Q^T-L} that 
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appear in Q are explicitly given by 

J --2 = + C-iC\ — diuir , 

g_2 = AX>Lia_f*I?o - 3AX>^^ia_ - C-ia+f*Cl 
+3£_ia+iasin0, 

T-L-2 = — AI?lj^Q;_f*£_i — 3AI?Lj^Qf_msin0 
-/:_ia+f*AX>li - 3£_ia+A, 

J--1 = lS.'DiD^_^+C\C-i — Qiujf, (7) 

g_^ = -Voa+f*M)l^-ZVoa+l^ + C\a_f*C_i 
+2iC\a-ia sin0, 

7^-1 = —I?oa+r*£2 + 327oC«+*asin0 — £2Q^_f*I?o 

with f = r + iacos0, a± = [3(f^M — rQ^) ± Q^r*] 
and we introduce the radial and angular Chandrasekhar 
operators, 

T>j = dr + + 2j-^-, Kr = am — (r^ + a^)oj; 

TTl 

Cj = de + Kg + j cot 9, Kg = — —- — aw sin d. (8) 

sin 9 

The complex conjugate of these operators, namely I?] 
and can be obtained from Vj and Cj via the replace¬ 
ment Kr —>■ —Kr and Kg —>■ —Kg, respectively. 
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